Chapter 17

Comparing Runge-Kutta methods.

In addition to the explicit Euler, implicit Euler, and trapezoidal methods introduced in the previous chapter, in this chapter we learned an additional class of explicit numerical methods for approximating solutions to first-order IVP’s, called Runge-Kutta methods.

Consider a (possibly non-autonomous) first-order IVP systems of the form \[ \frac{d\x}{dt} = \mathbf{f}(t,\x),\,\,\, \x(t_0) = \x_0. \] The goal is to approximate $\x(t)$ on the interval $t\in[t_0,t_0+T]$ for a chosen time duration $T > 0$. For this IVP, the modified Euler method (also called the midpoint method) is \[ \x_{k+1} = \x_k + h \mathbf{f} \left(t_k + \frac{1}{2}h, \x_k + \frac{1}{2}h \mathbf{f}(t_k, \x_k)\right) \] and the improved Euler method (also called Heun's method) is \[ \x_{k+1} = \x_k + \frac{1}{2}h \left( \mathbf{f}(t_k,\x_k) + \mathbf{f}(t_k+h, \x_k+h\mathbf{f}(t_k,\x_k)) \right). \]

For this IVP, the (fourth-order) Runge-Kutta method (also called "RK4") is \begin{equation}\label{rk4m} \x_{k+1} = \x_k + \frac{1}{6}h\mathbf{m}_1 + \frac{1}{3}h\mathbf{m}_2 + \frac{1}{3}h\mathbf{m}_3 +\frac{1}{6}h\mathbf{m}_4, \end{equation} where the $\mathbf{m}_j$'s are defined by the following recursive procedure: $$ \mathbf{m}_1 = \mathbf{f}(t_k,\x_k), \,\,\, \mathbf{m}_2 = \mathbf{f}\left(t_k+\frac{h}{2},\x_k+\frac{h}{2} \mathbf{m}_1\right),$$ $$ \mathbf{m}_3 = \mathbf{f}\left(t_k+\frac{h}{2},\x_k+\frac{h}{2}\mathbf{m}_2\right), \,\,\, \mathbf{m}_4 = \mathbf{f}(t_k+h,\x_k+h\mathbf{m}_3).$$

Below we give two examples of how the solutions of the three new methods (midpoint, Heun, and RK4) as well as the three old methods (explicit Euler, implicit Euler, and trapezoidal) approach to the exact solution. You can adjust the time step size using the slider. You can show or hide a curve by clicking on its legend. As you vary $h$, you will sometimes see blue exact solution curve move too: that curve really is not changing; what is happening is that the vertical axis scale sometimes changes to accommodate a smaller range of needed values for the output of the numerical methods. When you hover over the box, typically all that shows are black and blue number rectangles corresponding to coordinates of a point on the blue exact solution curve. When the black number is an integer multiple of $h$, then the output of the numerical methods also appears.

Comparing numerical methods with stiff ODE’s.

As discussed in this chapter, an ODE whose solutions converge very rapidly to a specific solution whatever the initial conditions (informally: solutions involve a rapidly decaying “transient” term, such as $C e^{-5x}$ for $y' = -5(y-\cos x)$) and an ODE system whose solutions involve multiple contributions that occur on time scales that differ by orders of magnitude is called a stiff ODE.

There is no precise mathematical definition of stiffness, but entire books have been written about numerical algorithms that handle stiffness with a reasonable balance between accuracy and efficiency. A common feature of stiff ODE's, beyond the presence of very different time scales within the same problem, is that “explicit” methods (such as the explicit Euler method) tend to perform very poorly. The reason that stiffness was overlooked prior to the advent of computers is that it is a phenomenon specific to working with actual step sizes $h > 0$ rather than a feature (such as stability) connected to behavior in the limit as $h \to 0$.

Same as in the textbook, here we illustrate through a study of the IVP \begin{equation}\label{basicstiff} \frac{dx}{dt} = -\lambda x,\,\,\, x(0) = 1 \,\,\,\mbox{ with } \,\,\, \lambda > 0 \end{equation} how to quantify the smallness of $h$ needed to overcome stiffness, and to show with examples that certain methods perform better than others. We will use a step size $h > 0$ and work on the time interval $[0, +\infty)$. Lest this ODE seem too special, we recall that through the method of eigenvectors we can generally regard most linear ODE systems as being a collection of many decoupled ODE's of the form $x'_j = -\lambda_j x_j$ (by the magical substitution). So phenomena that arise for this can be reasonably representative of what happens for typical large linear ODE systems (and even for large non-linear ODE systems near stationary points, through the process of linearization). In particular, the lessons we learn will be adapted to systems of linear ODE's by using eigenvalues of a derivative matrix in the role of $\lambda$.

Below we visualize this IVP. You can change the parameter $\lambda$ and adjust the step size $h$ using the sliders. You can show or hide a curve by clicking on its legend. As you vary $h$, you will sometimes see blue exact solution curve move too: that curve really is not changing; what is happening is that the vertical axis scale sometimes changes to accommodate a smaller range of needed values for the output of the numerical methods. As you vary $h$, observe that implicit methods are more resistant to stiffness issues (which appear when $h$ isn’t too small). When you hover over the box, typically all that shows are black and blue number rectangles corresponding to coordinates of a point on the blue exact solution curve. When the black number is an integer multiple of $h$, then the output of the numerical methods also appears.

Comparing different numerical methods for your favorite ODE or ODE system.

This widget will allow you to obtain numerical solutions to first-order initial value problems of your own choice, and this widget does the same thing for first-order 2-dimensional ODE systems.